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OO ■ Quantum simulation can beat current classical computers with minimally a 

(N ; 

. few tens of qubits and will likely become the first practical use of a quan- 

' 

^ '. tum computer. One promising application of quantum simulation is to at- 

^ ] '. tack challenging quantum chemistry problems. Here we report an experimen- 

^ '. tal demonstration that a small nuclear-magnetic-resonance (NMR) quantum 

■ computer is already able to simulate the dynamics of a prototype chemical 

a ■ 

reaction. The experimental results agree well with classical simulations. We 
conclude that the quantum simulation of chemical reaction dynamics not com- 
putable on current classical computers is feasible in the near future. 
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Introduction. In addition to offering general-purpose quantum algorithms with substantial 

speed-ups over classical algorithms (1) [e.g., Shor's quantum factorizing algorithm (2)], a quan- 
tum computer can be used to simulate specific quantum systems with high efficiency (i). This 
quantum simulation idea was first conceived by Feynman (4). Lloyd proved that with quantum 
computation architecture, the required resource for quantum simulation scales polynomially 
with the size of the simulated system (5), as compared with the exponential scaling on classical 
computers. During the past years several quantum simulation algorithms designed for individual 
problems were proposed (6-10) and a part of them have been realized using physical systems 
such as NMR (11-13) or trapped-ions (14). For quantum chemistry problems, Aspuru-Guzik et 
al. and Kassal et al. proposed quantum simulation algorithms to calculate stationary molecular 
properties (15) as well as chemical reaction rates (16), with the quantum simulation of the for- 
mer experimentally implemented on both NMR (17) and photonic quantum computers (18). hi 
this work we aim at the quantum simulation of the more challenging side of quantum chemistry 
problems - chemical reaction dynamics, presenting an experimental NMR implementation for 
the first time. 

Theoretical calculations of chemical reaction dynamics play an important role in under- 
standing reaction mechanisms and in guiding the control of chemical reactions (19,20). On 
classical computers the computational cost for propagating the Schrodinger equation increases 
exponentially with the system size. Indeed, standard methods in studies of chemical reaction 
dynamics so far have dealt with up to 9 degrees of freedom (DOF) (21). Some highly so- 
phisticated approaches, such as the multi-configurational time-dependent Hartree (MCTDH) 
method (22), can treat dozens of DOF but various approximations are necessary. So generally 
speaking, classical computers are unable to perform dynamical simulations for large molecules. 
For example, for a 10-DOF system and if only 8 grid points are needed for the coordinate rep- 
resentation of each DOF, classical computation will have to store and operate 8^° data points. 
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already a formidable task for current classical computers. By contrast, such a system size is 
manageable by a quantum computer with only 30 qubits. Furthermore, the whole set of data 
can be processed in parallel in quantum simulation. 

In this report we demonstrate that the quantum dynamics of a laser-driven hydrogen transfer 
model reaction can be captured by a small NMR quantum simulator. Given the limited number 
of qubits, the potential energy curve is modeled by 8 grid points. The continuous reactant- 
to-product transformation observed in our quantum simulator is in remarkable agreement with 
a classical computation based also upon an 8-dimensional Hilbert space. To our knowledge, 
this is the first explicit implementation of the quantum simulation of a chemical reaction pro- 
cess. Theoretical methods and general experimental techniques described in this work should 
motivate next-generation simulations of chemical reaction dynamics with a larger number of 
qubits. 

Theory. Previously we were able to simulate the ground-state energy of Hydrogen molecule 
{17). Here, to simulate chemical reaction dynamics, we consider a one-dimensional model 
of a laser-driven chemical reaction {23), namely, the isomerization reaction of nonsymmetric 
substituted malonaldehydes, depicted in (Fig. 1 A). The system Hamiltonian in the presence of 
an external laser field is given by 

H{t) = T + y + E{t) with E{t) = -iie{t). (1) 

In (Eq. 1), E{t) is the laser-molecule interaction Hamiltonian, ji = eqis the dipole moment 
operator, e{t) represents the driving electric field, T — /2m is the kinetic energy operator, 
and 

V=^{q- qo) + ~4^^^ g - qofiq + Qof (2) 

is a double-well potential of the system along the reaction coordinate. In (Eq. 2) is the barrier 
height, A gives the asymmetry of the two wells, and ±qo give the locations of the potential well 
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minima. See the figure caption of (Fig. IB) for more details of this model. 

We first employ the split-operator method {24, 25) to obtain the propagator U{t -\- dt, t) 
associated with the time interval from ttot + 5t. We then have 

U{t + St, t) « e-t^'5*/2g-t^(*+'5*/2)ft/2g-t™ 

X e-iE{t+St/2)St/2^-iVSt/2_ 

The unitary operator e"'™/'* in (Eq. 3) is diagonal in the momentum representation whereas 
all the other operators are unitary and diagonal in the coordinate representation. Such U{t + 
5t, t) can be simulated in a rather simple fashion if we work with both representations and 
make transformations between them by quantum Fourier transform (QFT) operations. To take 
snapshots of the dynamics we divide the reaction process into 25 small time steps, with 5t — 
1.5 fs and the total duration tf — 37.5 fs. The electric field of an ultrashort strong laser pulse is 
chosen as 

r £osin2(^); 0<t<si 
e{t) = <{ £o; si<t<S2 (4) 

with si = 5 fs and S2 — 32.5 fs. More details, including an error analysis of the split-operator 
technique, are given in the supplementary material. The reactant state at i = is assumed to 
be the ground-state |0o) of the bare Hamiltonian T + V, which is mainly localized in the left 
potential well. The wavefunction of the reacting system at later times is denoted by \'^{t)). The 
product state of the reaction is taken as the first excited state |0i) ofT -\-V, which is mainly 
localized in the right potential well. 

With the system Hamiltonian, the initial reactant state, the product state, and the propagation 
method outlined above, the next step is to encode the time-evolving wavefunction \i}}{t)) and the 
T, V, E{t) operators by n qubits. To that end we first obtain the expressions of these operators 
in representation of a set of N — 2^ discretized position basis states. The evolving state can 
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then be encoded as 



=mo(i) |0 • • • 00) + ... + m2n_i(i) |1 • • • 11) , 



(5) 



and as a result the system operators become 



T 



El 2 



.n 



(6) 



V 




(7) 



n 



Q 




(8) 



where (j = 1, 2, • • • ,n) is the PauU matrix and aj is the A'^-dimensional identity matrix. 
Because our current quantum computing platform can only offer a limited number of qubits and 
the focus of this work is on an implementation of the necessary gate operations under the above 
encoding, we have employed a rather aggressive 8-point discretization using n — 3 qubits. 
The associated diagonal forms of the T, V, and q matrices are given in the supplementary 
material. In particular, the end grid points are at g = ±0.8 A and the locations of other 6 grid 
points are shown in (Fig. IB). The eigenvalues of the ground and first excited states of the bare 
Hamiltonian treated in the 8-dimensional encoding Hilbert space are close to the exact answers. 
The associated eigenfunctions are somewhat deformed from exact calculations using, e.g., 64 
grid points. Nonetheless, their unbalanced probability distribution in the two potential wells 
is maintained. For example, the probability for the first excited state being found in the right 
potential well is about 80%. 

Experiment. In our experiment qubits 1,2, and 3 are realized by the ^^F, ^^C, and nuclear 
spins of Diethyl-fluoromalonate. The structure of Diethyl-fluoromalonate is shown in (Fig. 
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2A), where the three nuclei used as qubits are marked by oval. The internal Hamiltonian of this 

system is given by 



where Uj is the resonance frequency of the jih spin and Jjk is the scalar coupling strength be- 
tween spins j and k, with J12 — —194.4 Hz, /13 = 47.6 Hz, and 723 = 160.7 Hz. The relaxation 
time Ti and dephasing time T2 for each of the three nuclear spins are tabulated in (Fig. 2B). 
The experiment is conducted on a Bruker Avance 400 MHz spectrometer at room temperature. 

The experiment consists of three parts: (A) Initial state preparation. In this part we pre- 
pare the ground state |0o) of the bare Hamiltonian T + y as the reactant state; (B) Dynamical 
evolution, that is, the explicit implementation of the system evolution such that the continu- 
ous chemical reaction dynamics can be simulated; (C) Measurement. In this third part the 
probabilities of the reactant and product states associated with each of the 25 snapshots of the 
dynamical evolution are recorded. For the jth snapshot at tj = j5t, we measure the overlaps 



C(|Vfe)),l0o)) = |(0o|V'fe))l'andC(|V(i,)),|0i)) = through which the con- 



tinuous reactant-to-product transformation can be displayed. The main experimental details are 
as follows. Readers may again refer to the supplementary material for more technical explana- 



(A) Initial State Preparation. Starting from the thermal equilibrium state, firstly we create 
the pseudo-pure state (PPS) pooo = (1 — e)I/8 + e |000) (000| using the spatial average technique 
(26), where e fa 10~^ represents the polarization of the system and I is the 8 x 8 identity matrix. 
The initial state |0o) was prepared from pooo by applying one shaped radio-frequency (RF) 
pulse characterized by 1000 frequency segments and determined by the GRadient Ascent Pulse 
Engineering (GRAPE) algorithm (27-29). The preparation pulse thus obtained is shown in 
(Fig. 2C) with the pulse width chosen as 10 ms and a theoretical fidelity 0.995. Because the 
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(9) 



tions. 
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central resonance frequencies of the nuclear spins are different, (Fig. 2C) shows the RF field 
amplitudes vs time in three panels. To confirm the successful preparation of the state |0o)> we 
carry out a full state tomography and examine the fidelity between the target density matrix 
Po — 1 00 )( 00 1 and the experimental one Pexpi^)- Using the fidelity definition F(pi,p2) = 
Tr(pip2)/\/ (Tr(pf)Tr(pl), we obtain F[po, pexp{0)] = 0.950. Indeed, their real parts shown 
in (Fig. 4A) are seen to be in agreement. 

(B) Dynamical Evolution. The reaction process was divided into M — 25 discrete time 
intervals of the same duration 6t. Associated with the mth time interval, the unitary evolution 
operator is given by 

Um ~ Vst/2Est/2{tm)UQFTTstUQFTEst/2{tm)Vst/2, (10) 

where Uqft represents a QFT operation, and other operators are defined by Vst/2 = e~t^^, 
Tst = e-H™, and Est/2{tm) = ei^(*— i+'^*/2)^9t, with V, T, and q all in their diagonal repre- 
sentations. Such a loop of operations is m-dependent because the simulated system is subject 
to a time-dependent laser field. The numerical values of the diagonal operators Tst, V5t/2 and 
Est/2 are elaborated in the supplementary material. A circuit to realize Uqft and a computa- 
tional network to realize the Um operator are shown in (Fig. 3). 

Each individual operation in the Um loop can be implemented by a particular RF pulse se- 
quence applied to our system. However, in the experiment such a direct decomposition of Um 
requires a very long gate operation time and highly complicated RF pulse sequences. This 
bottom-up approach hence accumulates considerable experimental errors and also invites seri- 
ous decoherence effects. To circumvent this technical problem we find a better experimental 
approach, which further exploits the GRAPE technique to synthesize Um or their products with 
one single engineered RF pulse only. That is, the quantum evolution operator U {tj, 0), which is 
simulated by Y[m=i ^rn-, is implemented by one GRAPE coherent control pulse altogether, with 
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a preset fidelity and a typical pulse length ranging from 10 ms to 15 ms. For the 25 snapshots 
of the dynamics, totally 25 GRAPE pulses are worked out, with their fidelities always set to be 
larger than 0.99. As a result, the technical complexity of the experiment decreases dramatically 
but the fidelity is maintained at a high level. The task of finding a GRAPE pulse itself may be 
fulfilled via feedback learning control (19) that can exploit the quantum evolution of our NMR 
system itself. However, this quantum procedure is not essential or necessary in our experiment 
because here the GRAPE pulses on a 3-qubit system can be found rather easily. 

(C) Measurement. To take the snapshots of the reaction process at tj = j5t we need to 
measure the overlaps of C(\ip{tj)),\(f)o)) and C(|'0(ij)),|0i)). A full state tomography at tj will 
do, but this will produce much more information than needed. Indeed, assisted by a simple di- 
agonalization technique, sole population measurements already suffice to observe the reactant- 
to-product transformation. 

Specifically, in order to obtain C{\tlj{tj)) , |0o)) = Tr[p(ij)po] with p{tj) — 
we first find a transformation matrix R to diagonalize po, that is, pg = RpqR\ where pg is 
a diagonal density matrix. Letting p'{tj) — Rp{tj)R'^ and using the identity Tr[p(ij)po] = 
Tr[p'(ij)po], it becomes clear that only the diagonal elements or the populations of p'{tj) are 
required to measure Tr[p'(tj)pQ], namely, the overlap C{\ip{tj)) ,\4)q)). To obtain p'{tj) from 
p{tj), we simply add the extra R operation to the quantum gate network. The actual implemen- 
tation of the R operation can be again mingled with all other gate operations using one GRAPE 
pulse. A similar procedure is used to measure C(|V'(tj)),|0i))- 

The populations of p'{tj) can be measured by applying [7r/2]j^ pulses to the three qubits and 
then read the ensuing free induction decay signal. In our sample of natural abundance, only 
~ 1% of all the molecules contain a ^^C nuclear spin. The signals from the and ^^F nuclear 
spins are hence dominated by those molecules with the ^^C isotope. To overcome this problem 
we apply SWAP gates to transmit the information of the and ^^F channels to the ^^C channel 
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and then measure the qubit. 

To assess the difference between theory and experiment, we carry out one full state tomog- 
raphy for the final state density matrix att — tf. Because the GRAPE pulse is made to reach 
a fidelity larger than 0.995, the experimental density matrix Pexp{tf) is indeed very close to 
the theoretical density matrix ptheory{tf) obtained in an 8-dimensional Hilbert space, with a fi- 
delity F[pth^ory{^f)i Pexpi^f)] — 0.957. The experimental density matrix elements of the final 
state shown in (Fig. 4B) match the theoretical results to a high degree. With confidence in 
the experimental results on the full density matrix level, we can now examine the simulated 
reaction dynamics, reporting only the probabilities of the reactant and product states. (Fig. 4C) 
shows the time-dependence of the probabilities of both the reactant and product states obtained 
from our quantum simulator. It is seen that the product-to-reactant ratio increases continuously 
with time, with the probability of the product state reaching 77% at the end of the simulated 
reaction. At all times, the experimental observations of the reaction process are in impressive 
agreement with the smooth curves calculated theoretically on a classical computer. Further, the 
experimental results are also in qualitative agreement with the exact classical calculation using 
64 grid points (see Fig. IB). A prototype laser-driven reaction is thus successfully simulated 
by our 3 -qubit system. We emphasize that due to the use of GRAPE pulses in synthesizing the 
gate operations, our simulation experiment lasts about 30 ms only, which is much shorter than 
the spin decoherence time of our system. The slight difference between theory and experiment 
can be attributed to imperfect GRAPE pulses, as well as inhomogeneity in RF pulses and in the 
static magnetic field. 

Conclusion. Quantum simulation with only tens of qubits can already exceed the capacity 
of a classical computer. Before realizing general-purpose quantum algorithms that typically 
require thousands of qubits, a quantum simulator attacking problems not solvable on current 
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classical computers will be one conceivable milestone in the very near future. The realization 
of quantum simulations will tremendously change the way we explore quantum chemistry in 
both stationary and dynamical problems (15, 16). Our work reported here establishes the first 
experimental study of the quantum simulation of a prototype laser-driven chemical reaction. 
The feasibility of simulating chemical reaction processes on a rather small quantum computer 
is hence demonstrated. Our proof-of-principle experiment also realizes a promising map from 
laser-driven chemical reactions to the dynamics of interacting spin systems under shaped RF 
fields. This map itself is of significance because it bridges up two research subjects whose 
characteristic time scales differ by many orders of magnitude. 
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Fig. 1. Prototype chemical reaction and potential energy curve. (A) Isomerization reaction of 

nonsymmetric substituted malonaldehydes. (B) Upper panel: Potential energy curve, together 
with the eigenfunctions of the ground (red) and the first excited (blue) states. The main system 
parameters [(Eq. 2)] are taken from Ref. (25), with = 0.00625 E^, A = 0.000257 E^, and 
go = 1 oo- As a modification, the potential values for q approaching the left and right ends 
are increased sharply to ensure rapid decay of the wavefunction amplitudes, hi particular, this 
procedure increases the V value at g = ±0.8 A by a factor of 30. The six discrete squares 
shown on the potential curve and the two end points at g = ±0.8 A constitute the 8 grid points 
for our 3-qubit encoding. Lower panel: Numerically exact time-dependence of populations of 
the ground state (reactant state, denoted Pq) and the first excited state (product state, denoted 
Pi). 
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Fig. 2. (A) Molecular structure of Diethyl-fluoromalonate. The ^H, ^^C and ^^F nuclear spins 

marked by oval are used as three qubits. (B) System parameters and important time scales of 
Diethyl-fluoromalonate. Diagonal elements are the Larmor frequencies (Hz) and off-diagonal 
elements are scalar coupling strength (Hz) between two nuclear spins. Relaxation and dephas- 
ing time scales (second) Ti and T2 for each nuclear spin are listed on the right. (C) The GRAPE 
pulse that realizes the initial state |0o) from the PPS |000), with a pulse width 10 ms and a 
fidelity over 0.995. The (blue) solid line represents the pulse power in x-direction, and the 
(red) dotted line represents the pulse power in ^/-direction. The three panels from top to bottom 
represent the RF features at three central frequencies associated with the ^^F, ^^C and ^H spins, 
respectively. 



15 



Fig. 3. Upper panel: The network of quantum operations to simulate the chemical reaction 

dynamics, with the reactant state |0o)- The whole process is divided into 25 loops. The operators 
Tst, Vst and .^54/2 are assumed to be in their diagonal representations. Lower panel: H is the 
Hadamard gate and S, T are phase gates as specified on the right. Vertical lines ending with a 
solid dot represent controlled phase gates and the vertical line between two crosses represents a 
SWAP gate. 
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Fig. 4. Experimental tomography results and the reaction dynamics obtained both theoretically 

and experimentally. (A)-(B) Real part of the density matrix of the initial and final states of 
the simulated reaction. Upper panels show the theoretical results based on an 8-dimensional 
Hilbert space, and lower panels show the experimental results. (C) The measured probabilities 
of the reactant and product states to give 25 snapshots of the reaction dynamics. The (red) plus 
symbols represent measured results of C( |V'(ij)),|0o)) and the (blue) circles represent measured 
results of C(|V'(ij)),|0i)), both in agreement with the theoretical smooth curves. Results here 
also agree qualitatively with the numerically exact dynamics shown in (Fig. IB). 
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Supporting Online Material 

BACKGROUND ON QUANTUM DYNAMICS SIMULATION 

Let us start with the Schrodinger equation: 



Its formal solution can be written as 



iP{t)^U{t,to)ip{to), 
where the quantum propagator U{t, to) is a unitary operator and is given by 



U{t,to) = Texp 



(11) 



(12) 



(13) 



with T being the time ordering operator. There are a number of established numerical methods 
for propagating the Schrodinger equation, such as Feynman's path integral formalism (i), the 
split-operator method (2) and the Chebychev polynomial method (3, 4), etc. For our purpose 
here we adopt the split-operator method. 
The propagator U{t,to) satisfies 



U{t,to) = U{t,tN-l)U{tN-l,tN-2) 



U{h,to), 



(14) 



where, for example, the intermediate time points can be equally spaced with = rndt + to. 
For one of such a small time interval, e.g., from t^-i to t^ — t^-i + St, we have 



U{tm,tm-l) = Texp 



drHir) 



tm — 1 



exp 



h 



dTH{T) 



(15) 



where terms of the order {5tY or higher are neglected. For sufficiently small 5t, the integral in 
the above equation can be further carried out by a midpoint rule, leading to 



U{tm,tm-l) ^ exp 



--H{tm-i + St/2)5t 



(16) 
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This integration step has an error of the order of {5ty, which is still acceptable if the total 
evolution time is not large. Next we separate the total Hamiltonian into two parts: 

H{t) = Ho{t) + H'{t). (17) 

For example, Ho{t) is the kinetic energy part of the total Hamiltonian and H'(t) represents 
the potential energy part. In general Ho^t) and H'{t) do not commute with each other. The 
split-operator scheme (2) applied to (Eq. [T6l) then leads to 

U{tm,tm-l) ~ i^rn-l+5t/2)5tl2 ^-iHo(trn-l+5t/2)5t ^-IH' {trr,-l+St/2)&t/2 _ ^Jg-j 

The small error of this operator splitting step arises from the nonzero commutator between 
HQ{t) and H'{t), which is at least of the order of The advantage of the split-operator 

method is that each step represents a unitary evolution and each exponential in (Eq. [TSl) can 
take a diagonal form in either the position or the momentum representation. The error of this 
operator splitting step is in general smaller than that induced by the aforementioned midpoint 
rule integration in (Eq. [16]). Because in our work the total duration of the simulated chemical 
reaction is short, the above low-order approach already has a great performance. If long-time 
simulations with preferably larger time steps are needed in a quantum simulation, then one 
may use even higher-order split-operator techniques for explicitly time-dependent Hamiltonians 
{5,6). 

EXPERIMENTAL IMPLEMENTATION 

The experiment consists of three steps: (a) Initial state preparation, which is to prepare the 
ground state |0o) of the bare Hamiltonian T + V (representing the reactant state); (b) 25 dis- 
crete steps of dynamical evolution to simulate the actual continuous chemical reaction dynam- 
ics; (c) Measurement of the overlaps C{\ip{tj)) , |0o)) = \{<Po\i'itj))?' and C(|^/'(tj)) , |0i)) = 
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I {(f)i\tlj{tj)) p at tj — j5t, which is to show the transformation between the reactant and product 
states. 

A. Initial State Preparation 

To prepare the ground state |0o)> firstly we need to create a pseudo-pure state (PPS) from the 
thermal equilibrium state - a mixed state which is not yet ready for quantum computation pur- 

3 

poses. The thermal equilibrium state of our sample can be written as pther = J2 7i-^^' where 7^ 

i=l 

is the gyromagnetic ratio of the nuclear spins. Typically, 7c = 1, 7h = 4 and 7f = 3.7, with a 
constant factor ignored. We then use the spatial average technique (7) to prepare the PPS 

Pooo = ^I + e|000)(000|, (19) 

where e fa 10~^ represents the polarization of the system and I is the 8 x 8 unity matrix. The 
unity matrix has no influence on our experimental measurements and hence can be dropped. 
The pulse sequence to prepare the PPS from the thermal equilibrium state is shown in (Fig. 
SI A). In particular, the gradient pulses (represented by Gz) destroys the coherence induced by 
the rotating pulses and free evolutions. After obtaining the PPS, we apply one shaped pulse 
calculated by the GRadient Ascent Pulse Engineering (GRAPE) algorithm (8-10) to obtain the 
initial state |0o), with the pulse width 10 ms and a fidelity 0.995. In order to assess the accuracy 
of the experimental preparation of the initial state, a full state tomography (11) is implemented. 
The fidelity (12) between the target density matrix ptarget and the experimental density matrix 
Pexp is found to be 

F {ptarget, Pexp) = '^AptargetPexp) I \J (Tr (p?„^gej)Tr (p^^^) 

^ 0.95. (20) 

A detailed comparison between ptarget and Pexp is displayed in (Fig. SIB). 
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B. Dynamical Evolution 

To observe the continuous reactant-to-product transformation, we divide the whole time evolu- 
tion into 25 discrete steps. For convenience all variables here are expressed in terms of atomic 
units. For example, in atomic units e = 1, h = I, and 5t = 62.02. To exploit the split-operator 
scheme in (Eq. [T8l) . we let Hq = T, and H'{t) = V — eqe{t). The kinetic energy operator T 
is diagonal in the momentum representation, whereas the V operator and the dipole-field inter- 
action —eqe{t) operator are both diagonal in the position representation. We then obtain from 
(Eq.IH 

U{tm,tra-l) ^ VstE stUqFTTstUQFT^ (21) 

where the operators 

Vh = e"'^^, (22) 

2 

Tst = e-*™, (23) 
Eh = e'«^(*™-i+^*/2)f (24) 

2 

are assumed to be in their diagonal representations. 

To map U{tm,tm-i) to our 3-qubit NMR quantum computer, we discretize the potential 
energy curve using 8 grid points. Upon this discretization, operators Vst, Tst, and q become 
8x8 diagonal matrices. Numerically, their diagonal elements (denoted Vdiag, Tdiag, and qdiag. 
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respectively) are found to be 

Vdiag^ (293.78,-0.10,1.85,5.41, 

5.46,2.02,0.18,305.44) x lO'^; 
Tdiag = (0, 0.91, 3.63, 8.16, 

14.51,8.16,3.63,-0.91) x lO'^; 
Qdiag = (-1-51, -1.08, -0.65, -0.22, 

0.22,0.65,1.08,1.51). (25) 

To evaluate the Est operator, we also need to discretize the time-dependence of the electric 
field associated with the ultrashort laser pulse. For 25 snapshots of the reaction dynamics, we 
discretize the trapezoid-type electric field by 25 points, i.e., 

e{t) = [0.05, 0.42, 0.85, 1, ...1, 0.85, 0.42, 0.05] x 10"^ (26) 

The quantum gate network for the QFT operation that transforms the momentum repre- 
sentation to the coordinate representation is already shown in (Fig. 3). It consists of three 
Hardmard gates (H), three controUed-phase gates (S and T) and one SWAP gate (vertical line 
linking crosses). The Hardmard gate H is represented by the Hardmard matrix 

H^^A\ \V (27) 



V2 V 1 -1 

which maps the basis state |0) to :^(|0) + |1)) and |1) to :^(|0) — |1)). The phase gates S and 



T are given by 



and 



1 

i 



1 

e'^/^ 
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(28) 



(29) 



The matrix forni of the SWAP gate is 



SWAP = 



/ 1 \ 

10 

10 

\ 1 / 



(30) 



which exchanges the state of the F qubit and with that of the H qubit. 

GRAPE Pulses. Since there are hundreds of logical gates in the required network of quan- 
tum operations, a direct implementation of the gate operation network will need a large number 
of single-qubit rotations as well as many free evolutions during the single-qubit operations. 
This bottom-up approach will then accumulate the errors in every single-qubit operation. Con- 
siderable decoherence effects will also emerge during the long process. For example, we have 
attempted to directly decompose the network into a sequence of RF pulses, finding that the re- 
quired free evolution time for the 25 loops of evolution is more than 1 s, which is comparable 
to the T2 time of our system. To overcome these problems and to reach a high-fidelity quantum 
coherent control over the three interacting qubits, the unitary operators used in our experiment 
are realized by shaped quantum control pulses found by the GRadient Ascent Pulse Engineering 
(GRAPE) technique (8-10). To maximize the fidelity of the experimental propagator as com- 
pared with the ideal gate operations, we use a mean gate fidelity by averaging over a weighted 
distribution of RF field strengths to minimize the inhomogeneity effect of the RF pulses applied 
to the sample. 

For a known or desired unitary operator U tar get, the goal of the GRAPE algorithm is to find 
a shaped pulse within a given duration ttotai to maximize the fidelity 



where Ucai is the unitary operator actually realized by the shaped pulse and 2"- is the dimension 
of the Hilbert space. We discretize the evolution time ttotai into N segments of equal duration 



F=\Tr{Ul,^,tUoai)l'2 



(31) 
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— ttotai/N, such that Ucai — Un ■ ■ ■ U2U1, with the evolution operator associated with the 
7th time interval given by 

m 

Uj = exp{-iAt{nint + J2 ^k{j)nk)}. (32) 

k=l 

Here l-Lint is the three-qubit self-Hamiltonian in the absence of any control field, l-Lk represents 
the interaction Hamiltonians due to the applied RF field, and Uk{j) is the control vectors as- 
sociated with Hk- Specifically, in our experiment Mfe(j) are the time-dependent amplitudes of 
the RF field along the x and y directions, for the F-channel, the C-channel and the H-channel. 
With an initial guess for the pulse shape, we use the GRAPE algorithm to optimize Mfc(j) iter- 
atively until Ucai becomes very close to Utarget- More details can be found from Ref. {8). The 
GRAPE technique dramatically decreases the duration and complexity of our experiment and 
at the same time increases the quantum control fidelity. In our proof -of -principle demonstration 
of the feasibility of the quantum simulation of a chemical reaction, the task of searching for 
the GRAPE pulses is carried out on a classical computer in a rather straightforward manner. It 
is important to note that this technique can be scaled up for many-qubit systems, because the 
quantum evolution of the system itself can be exploited in finding the high-fidelity coherent 
control pulses. 

As an example, (Fig. S2) shows the details of one 15 ms GRAPE pulse to realize the 
quantum evolution from t = to = 15t (also combining the operations for initial state 
preparation and the extra operation R that is useful for the measurement stage). The shown 
GRAPE pulse is found by optimizing the frequency spectrum divided into 750 segments. The 
shown GRAPE pulse has a fidelity over 0.99. 

C. Measurement 

To simulate the process of a chemical reaction, it is necessary to measure the simulated reactant- 
to-product transformation at different times. To that end w e measure the overlaps of C{\ip{tm)) ■,\(t>o)) 
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and C{\il){tm)) , \4>i)) at tm — mSt. Here we first provide more explanations about how a diag- 
onalization method can reduce the measurement of the overlaps to population measurements. 
Without loss of generality, we consider the measurement C'dV'r) , |0o))- 

Cmtr)) , \M) = li'Pomh))]' = Tr[p(tv)po], (33) 

where /^(ty) = \^p{t-j) {'ip{t7)\ and po = |0o) (0o|- Let i? be a transformation matrix which 
diagonalizes po to a diagonal density matrix = RpoR} . Then 

Tr[p(t7)po] = Tr[Rp{t^)R^RpoR^] = Tr[p\tj)p'^i (34) 

where p'{ti) = Rp{t'j)RK Clearly then, only the diagonal terms (populations) of p'iti) are 
relevant when calculating Tr[p'(t7)pQ], namely, the overlap C{\ip{t-j) , |0o))- Hence only popu- 
lation measurement of the density matrix ^'(^7) is needed to obtain the overlap between \'4>{t7)) 
and the initial state. The GRAPE pulse that combines the operations for initial state preparation, 
for the quantum evolution, as well as for the extra operation R is shown in (Fig. S2). 

The three population-readout spectra after applying the GRAPE pulse are shown in (Fig. 
S3B-S3D), together with the ^^C spectrum for the PPS |000). The populations to be measured 
are converted to the observable coherent terms by applying a [7r/2]y pulse to each of the three 
qubits. For the ^^C spectrum shown in (Fig. S3C), four peaks from left to right are seen, with 
their respective integration results representing P(5) — P(7), P(6) — -P(8), P(l) — -P(3), and 
P(2) — P(4), where P{i) is the ith diagonal element of p^ti). Experimentally the four integrals 
associated with the four peaks in (Fig. S3C) are found to be -0.098, -0.482, -0.089 and 
-0.071, which are close to the theoretical values -0.047, -0.501, -0.114 and -0.041. Further 
using other readouts from the ^^F (see Fig. S3B) and ^H (see Fig. S3D) spectra as well as the 
normalization condition Yl^i=\ ^(0 — we obtain all the 8 populations and hence the overlap 
C{\'4>{t-j) , |0o))- The theoretical and experimental results for this overlap are 0.535 and 0.529, 
which are in agreement. A similar procedure is used to obtain C{\ip{tjn) , \4>i))- 
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The spectra of the and ^^F channel are obtained by first transmitting the signals of the ^^F 
and qubits to the ^"^C qubit using SWAP gates. With this procedure all the spectra shown in 
(Fig. S3) are exhibited on the ^^C channel. Indeed, because in our sample of natural abundance, 
only fa 1% of all the molecules contain a ^^C nuclear spin, the signals from the and ^^F 
nuclear spins without applying SWAP gates would be dominated by those molecules with the 
^^C isotope. 
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Fig. SI. Pulse sequence for the preparation of the PPS and a comparison between experimental 

and theoretical density matrix elements for the reactant state |0o). (A) Pulse sequence that 
implements the PPS preparation, with 9 — 0.647r and x(x, Y, Y) representing rotations around 
the x{—x, y, —y) direction. Gz represents a gradient pulse to destroy the coherence induced by 
the rotating pulses and free evolutions. and represent the free evolution of the system 
under Hint for 5.252 ms and 1.286 ms, respectively. (B) Comparison between measured density 
matrix elements of the initial state |0o) and the theoretial target density matrix elements based 
on the 8-point encoding. Both the real part and the imaginary part of the density matrix elements 
are shown. 
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Fig. S2. GRAPE pulse to simulate the quantum evolution of the reacting system from t — Oto 

t-j = 7St. The top, middle and bottom panels depict the time-dependence of the RF pulses ap- 
plied to the F-channel, C-channel and H-channel, respectively. The (blue) solid line represents 
the pulse power applied in the x-direction, and the (red) dotted line represents the pulse power 
applied in the ^/-direction. 
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Fig. S3. Measured spectra to extract the populations of the system density matrix before or after 

applying the GRAPE shown in (Fig. S2). (A) ^^C spectrum of the PPS |000) as a resuh of a 
[7r/2]j^ pulse applied to the ^^C qubit. The area (integral) of the absorption peak can be regarded 
as one benchmark in NMR realizations of quantum computation. (B)-(D) Signals from the 
^^F, ^'^C and qubits after applying the GRAPE pulse and a [7r/2]y pulse to each of the three 
qubits. All the spectra are exhibited on the ^^C channel through SWAP gates. The integration of 
each spectral peak gives the difference of two particular diagonal elements of the density matrix 

P'ih). 
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